#
# Download human lincRNA from ensembl/biomart
# 10/6/2011
# dario(dot)beraldi(at)roslin.ed.ac.uk
# See https://code.google.com/p/lincrna-pig/

library(biomaRt)
library(sqldf)
library(RODBC)

get.lincrna<- function(mart){
    ## Function to get lincRNAs
    lincrna<- getBM(	
                 attributes= c("ensembl_gene_id", "ensembl_transcript_id"),
                 filters= "biotype",
                 value= 'lincRNA',
                 mart= mart)
    return(lincrna)
    }

get.lincrna.seq<- function(mart, lincrna_id){
    ## Function to get cdna sequence of an ensembl transcript
    lincrna.seq<- getSequence(
            mart= mart,
            id= lincrna_id,
            type= "ensembl_transcript_id",     ## What type of id is it? Use listFilters() for valid names
            seqType= "cdna")
    return(lincrna.seq)
    }

## Write out as FASTA (if file exists, the sequence is appended!!)
lincrna.hs<- get.lincrna(useMart("ensembl", dataset= "hsapiens_gene_ensembl"))
lincrna.hs.seq<- get.lincrna.seq(useMart("ensembl", dataset= "hsapiens_gene_ensembl"), lincrna.hs$ensembl_transcript_id)
lincrna.mm<- get.lincrna(useMart("ensembl", dataset= "mmusculus_gene_ensembl"))
lincrna.mm.seq<- get.lincrna.seq(useMart("ensembl", dataset= "mmusculus_gene_ensembl"), lincran.mm$ensembl_transcript_id)

lincrna.file<- "F:/data/20110610_lincRNA/lincrna.fa"
exportFASTA(lincrna.mm.seq, file= lincrna.file)
exportFASTA(lincrna.hs.seq, file= lincrna.file)

## Execute blast of mouse and human lincRNA against pig genome:
blastoutfile<- "F:/data/20110610_lincRNA/lincrna.blastn"
sscrofa<- 'F:/data/ensembl/release-56/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa9.56.dna.toplevel.fa'
## system(paste('blastn ', '-query ', lincrna.file, '-db', sscrofa, '-task blastn', '-evalue 0.01', '-out', blastoutfile, '-outfmt 6'))
system(paste('blastn ', '-query ', lincrna.file, '-db', sscrofa, '-task blastn', '-evalue 1e-6 -best_hit_overhang 0.2 -best_hit_score_edge 0.05', '-out', blastoutfile, '-outfmt 6'))
#
# blastn -query lincrna.fa -task blastn -db ~/ensembl/release-56/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa9.56.dna.toplevel.fa -evalue 1e-6 -best_hit_overhang 0.2 -best_hit_score_edge 0.05 -out lincrna.blastn -outfmt 6
#

# -----------------------------------------------------------------------------
# After BLASTN...
# -----------------------------------------------------------------------------
conn<- odbcConnect(dsn= 'pgVitelleschi')
## See 20110617_lincrna.sql
no_hits<- sqlQuery(conn, "select case when qseqid like 'ENSMUS%' then 'mmusculus' ELSE 'hsapiens' END as species, count(*) as no_hits, qseqid from lincrna_blastn where evalue < 1e-9 group by qseqid order by no_hits desc;")
no_hits[1:10,]
hits_hist<- sqldf("SELECT *, CASE WHEN no_hits >= 50 THEN 50 ELSE no_hits END as hist FROM no_hits")

graphics.off()
windows(height= 10/2.54, width= 12/2.54)
par(las= 1, cex= 0.85)
hs<- hist(hits_hist$hist[no_hits$species == 'hsapiens'], breaks= 50)
mm<- hist(hits_hist$hist[no_hits$species == 'mmusculus'], breaks= 50)
plot(x= hs$mids[1:(length(hs$breaks)-1)], y= hs$counts, type= "o", lwd= 1, col= 'dodgerblue', pch= 19, xlab= '', ylab= '')
mtext(text= 'No. hits', side= 1, line= 2, cex= 0.85)
mtext(text= 'No. queries', side= 2, line= 3, las= 0, cex= 0.85)
points(x= mm$mids[1:(length(hs$breaks)-1)], y= mm$counts, type= "o", lwd= 1, col= 'firebrick4', pch= 19)
legend( "topright", legend= c('Human', 'Mouse'), col= c('dodgerblue', 'firebrick4'), lty= 'solid', pch= 19 )
title(main= 'Number of blast hits per query')
savePlot('M:/Documents/LabBook/LabBook_Figures/20110610_lincRNA_blasthits.emf', 'emf')

# -------------------------[ Tritume ]------------------------------------------

blastout<- read.table(blastoutfile, sep= '\t')
names(blastout)
sqldf("select count(distinct V1) from blastout;")
hist(blastout$V12, xlim= c(0, 100), breaks= 200)
writeClipboard(listAttributes(mart)[1][,]) ## Attribute code
writeClipboard(listAttributes(mart)[2][,]) ## 	 description
writeClipboard(listFilters(mart)[1][,])    ##
listFilters(mart)

sqldf("select count(distinct ensembl_gene_id) from lincrna") ## 1239 lincrna genes
sqldf("select count(distinct ensembl_transcript_id) from lincrna") ## 2932 lincrna transcripts

